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In the color glass condensate framework the saturation scale measured in deep inelastic scattering 
of high energy hadrons and nuclei can be determined from the correlator of Wilson lines in the hadron 
wavefunction. These same Wilson lines give the initial condition of the classical field computation 
of the initial gluon multiplicity and energy density in a heavy ion collision. In this paper the Wilson 
line correlator in both adjoint and fundamental representations is computed using exactly the same 
numerical procedure that has been used to calculate gluon production in a heavy ion collision. 
In particular the discretization of the longitudinal coordinate has a large numerical effect on the 
relation between the color charge density parameter g 2 /j, and the saturation scale Q B . Our result for 
this relation is Q B ~ 0.6g 2 /1, which results in the classical Yang-Mills value for the "gluon liberation 
coefficient" c ~ 1.1. 
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I. INTRODUCTION 

A useful description of the hadron or nucleus wave- 
function at high energy is to view the small x degrees of 
freedom as classical color fields radiated by classical static 
color sources formed by the large x degrees of freedom 
[2, 0, Hi ■ This description, known as the color glass con- 
densate (for reviews see e.g. [1, Q), provides a common 
framework for understanding both small x deep inelas- 
tic scattering (DIS) and the initial stages of relativistic 
heavy ion collisions, both of which can be understood 
in terms of Wilson lines of the classical color field. The 
cross section for small x DIS can be expressed in terms 
of the correlator of two Wilson lines in the fundamental 
representation (i.e. the dipole cross section), and the ini- 
tial condition for the classical fields that dominate the 
first fraction of a fermi of a heavy ion collision is deter- 
mined by these same Wilson lines. The inverse of the 
correlation length of these Wilson lines is known as the 
saturation scale Q s . The dipole cross section, can be 
determined from the dipole model fits to DIS data on 
protons [1, 0, H, [E EE] and nuclei [HI, E2 or extensions 
from the proton to the nucleus using a parametrization 
of the nuclear geometry [H El El, El EE E3- On 
the other hand there is a large body of both analyti- 
cal QJ, El HE Hi SI 111 [2l and numerical classical 

Yang-Mills (CYM) H M, H M HE HE EL M, Hi 

computations of the "Glasma" 34] fields in the initial 



stages of relativistic heavy ion collisions. 

The aim of this paper is to relate the parameters of 
these two types applications of color glass condensate 
ideas of the to each other more precisely. This is done 
by computing the Wilson line correlator in the McLerran- 
Venugopalan (MV) model [J, Hi] using exactly the same 
numerical method that has been used to compute the 
initial transverse energy and multiplicity in a heavy ion 
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collision. By doing this we can relate saturation scale 
Qs, whose numerical value can be determined from fits 
to DIS data, to the color charge density g 2 n that deter- 
mines the initial conditions for a heavy ion collision. The 
calculation relating these two parameters has been done 
analytically by several authors [H HE 113, H HE HE 52 • 
The procedure used to construct the MV model Wilson 
lines in this paper is the same as used in the numerical 
computations of the Glasma fields and differs from these 
analytical computations in two ways. Firstly, as noted 
also in Ref. [42|], the analytical computation is done by 
spreading out the color source in rapidity, while in the 
numerical computations this has not been done. We shall 
see that this introduces a factor of 2 difference in the ac- 
tual numerical relation between Q s and g 2 [i. Secondly 
the analytical result for the relation between Q s and g 2 [i 
depends logarithmically on an infrared cutoff that must 
be used in an intermediate stage of the computation, 
whereas in most of the numerical work the only such 
cutoff has been the size of the system. We shall also 
discuss the uncertainty arising from the non-Gaussian 
functional form of the Wilson line correlator and argue 
that it introduces an additional ambiguity at the 10% 
level. While the uncertainty from these aspects is para- 
metrically unimportant (a constant or a logarithm), they 
must still be better understood in order to increase the 
predictive power of the calculations. 

The logic of this paper is that, instead of treating the 
color charge density g 2 /i in the "Glasma" calculations as 
a free phcnomenological parameter, one should be able to 
relate it exactly, even the constant under the logarithm, 
to the saturation scale Q s measured in DIS experiments. 
When DIS measurements are used to determine the value 
of the saturation scale, choosing what treatment of the 
rapidity direction to use in the MV model is mostly a 
matter of convenience as long as the value of g 2 \x used 
is consistent with this chosen implementation. Let us 
note that our concern here is not as much the effects 
of high energy evolution on the wavefunction, but the 
parametrization of the region x ~ 0.01 relevant for cen- 



2 



tral rapidities at RHIC, which would be a reasonable ini- 
tial condition for solving the BK or JIMWLK equation. 
The glasma field configurations obtained are the boost 
invariant fields that serve as the background for study- 
ing things like instabilities in the classical [43|, [3, |45[ 
field and hi ghe r order contributions to particle produc- 
tion msi IS- 

We shall first introduce our notation for the Wilson 
line correlators in Sec. [ill Then our numerical results 
are presented in Sec. IIIII and their implications for the 
interpretation of some of the earlier phenomenological 
work discussed in Sec. HVl 



II. WILSON LINES AND GLASMA FIELDS 

Consider a high energy nucleus or a hadron moving 
along the £ + -axis. Its fast degrees of freedom can be 
considered as a classical color current 



J + = 9P{*t,x ), 



(1) 



which acts as a source to a classical color field represent- 
ing the slower partons 



(2) 



In the MV model the color charge density is taken to be a 
stochastic random variable with a Gaussian distribution. 
In covariant gauge Eq. ((2]) can be solved as 



A+(x ,x T ) = 



gp(x T ,x ) 



(3) 



The path ordered exponential of this field gives the Wil- 
son line in the fundamental representation 



-A+ 



(4) 



It is this quantity that will concern us in the following. 

The cross section for a virtual photon scattering off a 
high energy hadron or nucleus can be expressed in terms 
of the dipole cross section, which is determined by the 
correlator of two Wilson lines in the fundamental repre- 
sentation [BL [49l| 



C(x T -y T ) = (Tr[/t(x T )£/(y T )), 



(5) 



with the expectation value () evaluated with the distri- 
bution of the sources. 

The Wilson line in the adjoint representation is given 

by 

[/ Qb (x T ) = 2Tr [t a U\x T )t b U(x T )] • (6) 

The correlator of adjoint representation Wilson lines 

C(x T - y T ) = (U ab (xT)U ab (y T )} (7) 

is related to the gluon distribution of a nucleus [HI, H3, H3| 
(See Refs. (iTJ, [Hi], [Hj] for a discussion on the intricacies 




FIG. 1: The Wilson line correlator for g 2 fiL = 100 and 
different values of N y . Above: fundamental representation 
k 2 C(k), below: adjoint representation k 2 C(k). 



of defining a gluon distribution in this case.) With some 
algebra this the adjoint representation correlator can be 
related to a higher correlator of fundamental representa- 
tion Wilson lines 



C(x T - y T ) = (| Tr [t/ t (x T )[/(y T )] \ % ~ 1 



(8) 



which is the form we shall use to evaluate it numerically. 

The initial conditions for the glasma fields are deter- 
mined by the pure gauge fields (in light cone gauge) of 
the two colliding nuclei [l||, HTJ] . In terms of the Wilson 
line HI) the pure gauge field of one nucleus is 



A<(x T ) = -f/(x T )9,C/t(x T ), 
9 



(9) 



and the initial conditions for the glasma fields are given 
by the sum and commutator of the pure gauge fields 
of the two nuclei. In the numerical computation of the 
glasma fields [U H, M, H, H, M, [M| there has been 
no longitudinal structure in the source, and the Wilson 
lines have been constructed simply as 



C/(xt) = exp 



■ 9P{~xt) 



(10) 



with the transverse charge densities depending on a single 
parameter fi, independent of x~: 



(p a (x T )p b {y T ))=5 ab 5 2 {x T 



\ 2 2 

yrjg n ■ 



(11) 
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The analytical calculation [3J, |3J, |37J, |3J, IM |40j, |4l| 
of the Wilson line correlator requires that, unlike the nu- 
merical procedure in ^lH^HIiiSISH, the source 
be extended in the x~ direction 



(p a (x T ,x-)p b (y T ,y-)) = 

g 2 5 ab 5 2 (x T -y T )5(x- 



y-)fi 2 (x-). (12) 



With this longitudinal structure the Wilson line corre- 
lators can be computed analytically up to a logarithmic 
infrared cutoff that must be introduced in solving the 
Poisson equation ([3]). The result is 



C(x T ) 
C(x r ) 



^pg-g^XXr ln(m|x T |) 



(13) 



with 



X = g 4 / da; i?(x ). 



(14) 



The dimensions and Casimirs of the two representations 
in Eq. ^ are d A = N c 2 - 1, dp = N c , C A = N c and 
C F = (N c 2 - l)/2N c . It could be argued that the cutoff 
m should be ~ Aqcd- In any case, running coupling and 
confinement effects are not included in this treatment 
and the cutoff cannot be consistently defined within this 
calculation. When looking at length scales |xy| <C 1/m 
results depend very weakly on this cutoff; in the lattice 
calculation it can be replaced by the finite size of the 
lattice. It would be very tempting to identify /i 2 , the 
source strength of the delta function source, appearing in 
Eq. (fTT]) . with the integral over the spread distribution 
H 2 (x~) of Eq. (|14p . but as we will see in the following, 
this correspondence is not exact. 1 

Note that the form (fT3")) is compatible with the expec- 
tation that in the large N c limit the four point function 
in Eq. ((8]) factorizes into a product of two point functions 
and 



A', 



lim C(x T ) = C 2 (x T ) 



(15) 



III. NUMERICAL PROCEDURE AND RESULTS 

The Wilson lines used in the numerical calculation of 
the Glasma fields [H, H, S3, HE US M, U S3 are 
SU(3) matrices defined on the sites of a 2 dimensional 
discrete lattice corresponding to the transverse plane. As 
in most of these calculations, we shall consider a square 



It 



relatively easy to 



that the identification of 



J di~fi 2 (xr,n~) with /i 2 (xt) of Eq. {TTJ would be exact in 
the Abelian case or in the large N c limit in which the terms 
resulting from the noncommutative nature of p are suppressed. 



0.9 
0.8 
-£? 0.7 

a 

0.6 
0.5, 



1 1 1 1 1 


$ 




$ 


t 














* 






o 


N = 

y 


1, fund 










w — 

IN — 

y 


i, aaj 




- 







N = 

y 


10, fund 


- 






♦ 


N = 

y 


10, adj 






CD 


cc 


CD 






















0.2 0.4 



0.6 



0.8 



1.2 



FIG. 2: The lattice spacing dependence of the saturation 
scales Q s and \J~%^;Qs for g 2 [iL — 100 and different N y . The 
continuum limit is the g 2 fia — axis on the left. 
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FIG. 3: The dependence of the adjoint and fundamental 
representation saturation scales Q B and J -§-^Qs on g 2 ^L for 
g 2 /ia = 0.5 and different N y . 



lattice with periodic boundary conditions and an average 
color charge density g 2 ^ 2 that is constant throughout the 
plane. The Wilson lines are constructed as follows: On 
each lattice site x-r one constructs random color charges 
with a local Gaussian distribution 



( P £(x T ) ft b (y T )) = 5 ab 6 kl 6 2 (x T - y T ) 



(16) 



with the indices k,l = l,...,N y representing a dis- 
cretized longitudinal coordinate. The numerical calcu- 
lations so far have been done using N y — 1, whereas the 
derivation of the analytical expression of the correlator, 
Eq. (|13p are derived with an extended source, correspond- 
ing to the limit N y — > oo. Our normalization is chosen 
so that 

]T (p£(x T )pf(y T )) = S ab S 2 (x T - y T ) 5 V- (17) 

k,l 



4 



1.2 







o 


m 


= 0, fund 


1.1 




• 


in 


= 0, adj 


1.0 




A 


m 


2 

= 0.1 g H, fund 






A 


m 


= 0.1gV, adj 



0.9 

a 0.8 



8 A 



ft A 



J *1 1 



A A A * 



10 



N 



100 



FIG. 4: Dependence on N y of the saturation scales Q a and 
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The Wilson lines are then constructed from the sources 
(|16|) by solving a Poisson equation and exponentiating: 



N v 

n ° x p 

fe=i 



Pfe(x T ) 



(18) 



Here we have introduced an infrared regulator to for in- 
verting the Laplace operator. This is the same regulator 
as the one appearing in the analytical expression Eq. (|13[) . 
For large N y the charge densities p^in Eq. p^|) become 



small, and the individual elements in the product (|18[) 
approach identity. This is precisely the procedure that 
leads in the N y — > oo limit to the continuum path ordered 
exponential ([4]). 

To summarize, our calculation depends on the follow- 
ing parameters: 

• <? 2 p, determining the color charge density. 

• N yi the number of points in the discretization of 
the longitudinal (x~~ or rapidity) direction. 

• The infrared regulator to. When to = 0, as in 
most of the results presented, the Poisson equa- 
tion is solved by leaving out the zero transverse 
momentum mode. This procedure corresponds to 
an infrared cutoff given by the size of the system. 

• The lattice spacing a. 

• The number of transverse lattice sites N±, giving 
the size of the lattice L = N±a. 

Of the parameters a, g 2 fi and to, only the dimension- 
less combinations g 1 ! /xa and ma appear in the numerical 
calculation, and the continuum limit a — > is taken by 
letting N± — > oo so that g 2 fia — > and g 2 /iL = g 2 /iaN± 
remains constant. What we are looking at is relatively in- 
frared quantity and thus should converge very well in the 



continuum limit. Based on the analytical calculation we 
may expect a logarithmic dependence of the saturation 
scale on g 2 fi/m or, for to = 0, on g 2 [iL. 

By Fourier transforming the Wilson lines we can than 
construct the momentum space correlators in the adjoint 
and fundamental representations, C(\z.t) and C(k*r) re- 
spectively. These correlators, averaged over the polar 
angle, for different values of N y are plotted in Fig. [T] as 
a function of 



^sin 2 (ha/2). 



(19) 



For small momenta the correlators look like Gaussians, 
which is the form used in the "GBW" fit of DIS data 
m Refs. [1,0,1- For large momenta there is a power 
law tail 1/kr 4 that differs from the original GBW fits, 
but resembles more closely the form required to match 
smoothly to DGLAP evolution for large Q 2 [53| ■ 

We define the numerically measured saturation scales 
as follows. The scale Q s is determined by the adjoint 
representation Wilson line correlator as the momentum 
fcmax corresponding to the maximum of /c 2 C(k*r). This 
normalization in terms of the adjoint representation cor- 
responds to that of Refs. [13, Hy|. Similarly, from the 
maximum of the fundamental representation correlator 
k 2 C (kr) we define the fundamental representation sat- 
uration scale Qs as Q 2 = fc^. Our definition of the 
saturation scale is not sensitive to the exact shape of the 
correlator for very large or small transverse momenta, 
and for a Gaussian correlator it reproduces the GBW 
saturation scale as 1/Rq = Q 2 - The saturation scale is 
expected to scale according to the Casimir of the repre- 
sentation, meaning Q 2 w §jQs 2 - In the plots (Figs. [2 02 

2] and O we shall rescale Qs by this color factor to make 
the validity of this scaling clearer. 

We first check the lattice spacing dependence of our re- 
sult. Figure shows that, as expected, the ratio Q s /g 2 [i 
depends in fact so little on the lattice spacing that we 
will in the following not perform any continuum extrap- 
olation for this quantity. The dependence of Q s on the 
lattice size through the combination g 2 ixL (without the 
additional infrared cutoff m) is shown in Fig. [3] The 
values used in the numerical c omp utations of the glasma 
fields [H, [H, H3, [H, IH M, M, HI H| correspond to 
N y = 1 and g 2 /j,L ~ 100 in Fig. d with Q s « 0.57#V 

Figure [4] shows the dependence of g 2 fj,/Q s on the num- 
ber of points used to discretize the longitudinal direction, 
N y . When m = 0, i.e. the infrared singularity is regu- 
lated only by leaving out the zero mode, there is approx- 
imately a factor of two difference between Q s = 0.57<? 2 /i 
for N y = 1 (the numerical CYM prescription) and Q s w 
1.15g 2 /i for N y — > oo (the analytical computation of the 
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FIG. 5: The same fundamental representation Wilson line 
correlator as in Fig. [1] plotted as a function of the scaling 
variable k/Q s . Above: fundamental representation k 2 C(k) 
vs. k/Qs, below: adjoint representation k 2 C(k) vs. k/Qs- 
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dipole cross section) 2 . When a regulator m = 0.1<? 2 /z is 
introduced the dependence on N y is weaker, which can 
also be seen in Fig. [Sj In Fig. [5] we show the same corre- 
lators as in Fig. [T] as a function of k/Q s instead of k/g 2 /j,. 
One sees that the correlator has a scaling form indepen- 
dent of N y ; from which only the N y — 1 result deviates 
slightly. This suggests that, as argued in Sec. U once 
the appropriate relation between Q s and g 2 /! is used, the 
physical results depend very little on N y . Thus no sig- 
nificant change to the numerical CYM results should be 
expected if the calculations were repeated using a differ- 
ent treatment of the longitudinal coordinate in the source 
P- 

Explicitly regulating the infrared behavior with a mass 
scale m makes it possible to compare the numerical re- 
sult to the analytical one of Eq. (flUl) . If one introduces 
an infrared scale m as in Eq. (JTDJ) and replaces In (m |x<r|) 
with — ln(g 2 /i/m) in the coordinate space correlator it 
becomes a Gaussian. Fourier transforming this one ob- 



2 Because the initial energy density e of the glasma is proportional 
to Qs 4 /g 2 , this factor of 2 could be an explanation of the factor 
of 16 difference in e/(g 2 fi) 4 observed in Ref. O . 



FIG. 7: Dependence on the regulator m of the coordinate 
space saturation scale Q s , coord, for g 2 fiL — 100. The solid 
line is the expected logarithmic dependence, Eq. (|21[) with 
the constant C = 0.616 and the dashed one with C — 0. 
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Because of the replacement \x.t\ ~ l/sV there is still 
an uncertainty in the constant term. In Fig. [S] we plot 
the numerical result for Q s /g 2 /j, as a function of m/g 2 fi 
compared to the estimate 



Ca 
2tt 



(21) 



with values C = | 



In 2 - 7e rj 0.616 and C = 0. In 
an intermediate range of m) g 2 [x and for a large enough 
value of N y (recall that the analytical result corresponds 
to N, 

normalization different 

Another common way to define the saturation scale 
is in terms of the coordinate space correlator C(xt), 
because this is the object appearing in the calculation 



y —r oo) the behavior of Q s /g A* is similar, but the 
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FIG. 8: The ratio of the two definitions of the saturation 
scale, the coordinate space definition (Q s , coord.) and the mo- 
mentum space definition used in most of this paper (Qs)- 



of most observables in DIS. Kowalski and Teaney [ldj 
define the saturation scale Q s . coord, from the condition 
that C(x T ) = d A e' 1/2 at x 2 . = 2/Q 2 coord . Note 
that the definition in Ref. [TtJ where the same IPsat 
model is used differs slightly: C(xt) = d A e~ x ^ at 



Xy = 1/Q 2 coor( j • This definition can also be used in the 
numerical CYM computation, most straightforwardly by 
Hankel-transforming the correlator C(|kr|) back into co- 
ordinate space. As shown in Fig. using this definition 
is closer to the analytical estimate Eq. (l2"Tj) . The ratio 
of Qs, coord, to our original definition of Q s for different 
values of g 2 [iL and N y is plotted in Fig. [5J The differ- 
ence between the two definitions is of the order of 10% 
with small variations. One must emphasize here that for 
an exactly Gaussian Wilson line correlator (the GBW 
form) the two definitions would be equal. They differ in 
the MV model, because the correlator is not Gaussian. 
Thus if one tries to determine g 2 'fx from a comparison to 
the experimental DIS data using GBW-type fits, which 
is one of the alternatives we consider in the next sec- 
tion, the ambiguity in the definition of Qs leads to a 10% 
uncertainty in the value of g 2 [i- 



IV. DISCUSSION 

Let us finally use the results of the previous section for 
Q s /g 2 H and studies of DIS data to estimate the relevant 
value of g 2 fi for RHIC physics. In deep inelastic scatter- 
ing the variables x and Q 2 are precisely defined, and the 
saturation scale is a function of x, typically Q 2 ~ x~ x 
with A ~ 0.3. In the context of a heavy ion collision one 
is in fact, at a fixed energy and rapidity, summing up glu- 
ons produced at different transverse momenta and thus 
related to partons of different x in the nuclear wavefrmc- 
tion. The value of x at which to evaluate the satura- 
tion scale must therefore be some kind of effective x e f{, 





Qs H.hic 




Naive A 1/3 


1.7 GeV 


3.0 GeV 


A 1/A [11] 


1.4 GeV 


2.5 GeV 


C24 4/9 [12] 


1.3 GeV 


2.2 GeV 


IPsat, ~ CA 1/3 In A [17] 


1.1 GeV 


2.0 GeV 



TABLE I: Results for the adjoint representation saturation 
scale from extrapolations of DIS data to RHIC central rapid- 
ity kinematics. 



depending on the typical transverse momentum of the 
produced gluons, 



XeS ' 



(pt) ^ Qs_ 



(22) 



This introduces an additional uncertainty into our at- 
tempt to determine the color charge density based on 
the deep inelastic scattering data; by varying \Q s /\fs < 
x c ff < 2Q s /^/s we get an uncertainty of the order of 5%. 

A simple starting point for our estimate is the GBW 
fit @, 0] , where the proton saturation scale (in the fun- 
damental representation, which is convenient for DIS) is 
parametrized as 



l 



p2 

GBW 



foi.x/x ) 



(23) 



The result of the fit including charm quarks gives A 
0.277 and Xq = 0.41 • 10~ 4 , with the one redundant pa- 
rameter chosen as Qo = 1 GeV, while the fit without 
charm makes Q s approximately 30% larger. 

This result must then be extended to finite nuclei. 
Let us denote the nuclear modification of Q s by g(A) = 
QsA_/Qsp- The most straightforward theoretical expec- 
tation for the nuclear dependence would be g(A) — 
A 1 / 3 . Freund et al. [Hj perform a fit of the form 
g(A) = A s to the available nuclear DIS data and ob- 
tain (5 = 1/4. Taking into account modifications to the 
A 1 / 3 behavior of the nuclear radius leads Armesto et 
al. [12| to consider a fit of the form g(A) ~ AR 2 /R\ = 

C [A/ (A 1 / 3 - 0.77 A' 1 / 3 ) 2 ] 36 with the result S « 0.42 
and C ~ 0.5. Although for asymptotically large nu- 
clei this would imply g(A) ~ 7l 42 , for the physical case 
A < 200 the nuclear modification factor g(A) obtained 
in Ref. [l2j is actually less than 2L 1 / 4 . 

A more detailed description of the saturation scale in a 
nucleus can be obtained by the IPsat model [l(| EH, E3] • 
HERA data and the DGLAP equations are used to 
parametrize the dipole cross section for a proton. Tak- 
ing into account the fluctuations in the positions of the 
nucleons in the nucleus within a realistic nuclear geom- 
etry leads to a nuclear dipole cross section, from which 
also the saturation scale can be determined. As shown 
in Ref. [l7| this picture leads to a good parameter free 
description of all the existing small x eA data. The re- 
sult is a more realistic picture of an impact parameter 



7 



dependent saturation scale also influenced by DGLAP 
evolution, where the nuclear geometry leads to a g(A) 
that can roughly be understood as a C A 1 / 3 -like depen- 
dence (with C < 1) enhanced by a logarithmic increase in 
A resulting from the DGLAP evolution. Because scatter- 
ing off nuclei is less dominated by the dilute edge than in 
the proton, the typical Q s (conveniently taken as corre- 
sponding to 6 mo( j. , the median impact parameter in deep 
inelastic scattering) is closer to the maximal Q a in the 
nucleus than in the proton. We shall use here for gold 
the value at 6 me d. ~ 4.2 fm and at at x = 0.005 » Q s /\/~s- 

Table U summarizes the estimated saturation scales for 
calculating the classical field at central rapidity in RH1C 
based on the different fits explained above. The table 
also shows the corresponding values of g 2 /j, — Q s /0.57. 

To set this result in perspective let us briefly recall the 
result of the numerical CYM computations. The energy 
and multiplicity per unit rapidity can be parametrized as 



diV 
<h) 
dE T 
drj 
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(g 2 ^R 2 A 
g 2 



N 



Je 



(24) 
(25) 



where the numerical result is fp « 0.25 and /jy ~ 0.3 
(see in particular Refs. [1^, [3f[ for the result). This 
leads to the estimate that, for central rapidity at RHIC, 
1.3 GeV < g 2 n < 2.1 GeV. The lower limit comes 
from the requirement that the energy per unit rapidity 
has to be at least as large initially as it is finally. In 
boost invariant hydrodynamical ("Bjorken") flow energy 
is transferred from central rapidities to the fragmenta- 
tion region by p dV work thus decreasing the energy at 
central rapidities, decreasing dE?/ dy. It is very hard 
to imagine a process that would increase the energy at 
midrapidity and thus the measured final transverse en- 
ergy gives a lower limit to the initial energy and to g 2 \x. 
The upper limit follows from the requirement that the 
number of gluons in the initial state should be less or 
equal to that of hadrons in the final state. In ideal hy- 
drodynamical flow the two are related by entropy con- 
servation, and nonequilibrium processes should increase 
entropy and consequently the multiplicity during the evo- 
lution, not decrease it. The measured hadron multiplic- 
ity thus gives an upper limit on the initial multiplicity 
and g 2 fJ-- Quark pair production [55| or in general 
higher order processes would generically increase the ini- 
tial multiplicity for a given g 2 [i and thus decrease the 
upper limit for g 2 /j, below 2 GeV. The only overlap re- 
gion between these estimates and the DIS based ones in 
Table Q] is around g 2 /i w 2 GeV. This is also close to the 
estimate, based on large x parton distribution functions, 
of Re f . 12011 that was used in the CYM calculations of 
Refs. Halllii- On the DIS side the value g 2 [i w 2 GeV 
agrees very well with the estimate using the IPsat model 
(see table IT]), but our present understanding of gluon pro- 
duction seems to be in contradiction with the other, more 
naive, DIS fits. On the ion-ion collision side, to a large 
initial gluon multiplicity and rapid equilibration, leaving 



little room for higher order effects or additional gluon 
production during thermalization. 

One feature of the numerical CYM calculations has 
been the apparent small fraction of the gluons in the ini- 
tial nuclear wavefunction that are "freed" in the collision 
[50L [5^. The liberation coefficient c, introduced by A. 
Mueller [57| . is defined by writing the produced gluon 
multiplicity as 



dN 



— c 

d 2 x.T dy 2~K 2 a s 



With Eq. |2U this leads to 



2C f V Q 



g 2 v 



(26) 



(27) 



The original expectation was that c should be of order 
unity. The analytical calculation by Y. Kovchegov (50j 
gave the estimate c « 2 In 2 s» 1.4. Using the formula 
(|2T)|) for the ratio Q s /g 2 /i led to the interpretation [56|, 
l58j that the CYM result would be c « 0.5. We now 
see that when Q s /g 2 /j, is computed consistently with the 
numerical calculation the resulting CYM value for the 
liberation coefficient is c rs 1.1. We must emphasize that, 
because c is defined in terms of the physical multiplicity 
and the physical correlation length Q s , there is no large 
logarithmic or N y uncertainty in the result c ~ 1.1. The 
non-Gaussianity of the MV model correlator, as seen in 
the differing coordinate and momentum space results for 
Q s , does introduce an ambiguity at the 10% level. 

Let us summarize the major sources of error in estimat- 
ing the relevant value of the saturation scale for RHIC 
physics from the DIS data. We have already mentioned 
the questions of the Wilson line correlator not being ex- 
actly of the GBW form, the exact value of x to use and 
the considerable variance in the estimates of A depen- 
dence of Q s . It is also possible that including a more real- 
istic description of the transverse coordinate dependence 
of the saturation scale 53, E3, HO, HI] in the CYM 
calculation will have an impact on the gluon multiplicity 
and energy in an ion-ion collision, modifying our previous 
discussion. The solution to the problems related to the 
shape of the correlator can be solved by using the actual 
solution of the BK or JIMWLK equations to understand 
both DIS data (as is done in Ref. [§]) and to calculate 
the Glasma fields. Confirming the calculations like that 
of Ref. [17[ relating the saturation descriptions of the pro- 
ton and a nucleus will require more experimental input in 
the form of more data on small x DIS on nuclei. Finally 
and perhaps most importantly, the influence of quantum 
corrections and instabilities of small rapidity-dependent 
fluctuations is not yet understood quantitatively. 
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